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Application of three body stability to globular clusters I: 
the tidal radius 
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ABSTRACT 

The tidal radius is commonly analytically determined by equating the tidal field 
of the galaxy to the gravitational potential of the cluster. This gives a radius at which 
stars can move from orbiting the cluster centre to independently orbiting the galaxy. 
In this paper the tidal radius of a globular cluster is estimated using a novel approach 
from the theoretical standpoint of the general three-body problem. This is achieved 
by an analytical formula for the transition radius between stable and unstable orbits 
in a globular cluster. 

A stability analysis, outlined by iMardlinel (|2008l h is used here to predict the 
occurrence of unstable stellar orbits in the outermost region of a globular cluster 
in a distant orbit around a galaxy. It is found that the eccentricity of the cluster- 
galaxy orbit has a far more significant effect on the tidal radius of globular clusters 
than previous theoretical results have found. A simple analytical formula is given for 
determining the transition between stable and unstable orbits, which is analogous to 
the tidal radius for a globular cluster. This tidal radius estimate is interior to previous 
estimates and gives the innermost region from which stars can random walk to their 
eventual escape from the cluster. 
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1 INTRODUCTION 

The tidal radius of a globular cluster (GC) is defined as 
the point at which stars will escape the cluster's potential 
well and become part of the galactic halo. It is typically 
calculated by considering the equipotential poi nt between 
the cluster and t he tidal field of the galaxy fsee lKinglll962l ; 
iRead et al.ll2006T l. The aim of this paper is to determine the 
boundary between tidally stable and unstable orbits in a 
cluster potential, where an unstable orbit refers to a star 
orbiting inside the globular cluster that will eventually es- 
cape the cluster. This stability boundary is analogous to the 
tidal radius of a glo bular cluster and is predicted using the 
stability analysis of Mardling |200cf ) for GCs on eccentric 
galactic orbits and with arbitrary mass ratios. 

For the purposes of applying this stability analysis, the 
star-cluster-galaxy system is approximated as follows. A star 
of mass rrn orbits a particle representing the total cluster 
mass Mc which itself orbits the galaxy, taken as a particle 
of mass Mq- Each of these particles is treated as a point 
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mass so that we can treat the system as the general three- 
body problem. 

Simplifying the star-cluster-galaxy system to a single 
three-body problem ignores the effect of mutual interactions 
between stars in the cluster, which results in two-body relax- 
ation being ignored. This means that in a real cluster stars 
will be able to diffuse over the predicted tidal radius and 
escape the cluster from orbits that were initially in stable 
regions. However the very long relaxation timescale in these 
regions means that this effect is negligible. 

By approximating the cluster potential as a point mass 
in the stability analysis we are assuming that stars spend 
most of their time in the outermost regions of the cluster. 
Similarly by treating the galaxy as a point mass the discus- 
sion is limited to GCs, which spend most of their time at 
distances greater than approximately 5 kpc. The assumption 
of a point mass potential for the galaxy for distant cluster 
orbits is also justified since the mass inside 6.4 kpc is con- 
sistent with that of a point mass of roughly 1O 11 M0, based 
on the reasonably well kno wn orbits of the globular clusters 
NGC 2419 and NGC 7006 (|Bellazzinill2004T ). 

This paper is structured in the following way. A brief 
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overview of the Mardling stability criterion is presented in 
Section [2] This stability analysis, applied to the three-body 
system corresponding to a point mass approximation for the 
star, cluster and galaxy system, is presented in Section [3] A 
model cluster is simulated in isolation to get the distribution 
of star-cluster orbital eccentricities in Section 2] Section [S] 
applies the stability boundary to calculate the tidal radius 
for a case study globular cluster of mass M c = 10 6 A'/ Q , which 
is compared to previous theoretical work from the literature 
in Section [6] Conclusions and a discussion of the observa- 
tional and numerical simulation implications of this work 
are summarised in Section [7] 



2 MARDLING STABILITY CRITERION 

We first summarise the Mardling stability criterion (MSC) 
and outline the algorithm used to determine the stability of 
any three - body system. The criterion for the MSC is given in 
iMardlinsl (120081 ). and interested readers are referred to this 
study for details which have been omitted here. The criterion 
below covers the more general case where th e inner and oute r 
orbits are relatively inclined by an ang le I (|Mardlindl20"llh . 
The inner orbit refers to the orbit of the binary composed 
of point masses mi and 7712, while the outer orbit refers to 
the orbit of 7713 with the centre of mass associated with the 
combined mass m\i = mi + mi. 

We will use the common terminology of a system be- 
ing in resonance if the ratio of the outer to inner periods 
(a — To/Ti) is within a distance Act of a particular (n:n') 
resonance. When a system has initial conditions such that it 
resides in two resonances simultaneously then it is possible 
for the resonance angles to librate inside either resonance at 
any time. This jumping between resonances can occur at any 
time and means that two three-body systems with initially 
very close conditions will quickly diverge once one jumps to 
a different resonance, leaving the other behind. This sensi- 
tivity to changes in the initial conditions is characteristic of 
chaotic motion and is used by the MSC to predict unstable 
systems. 

For a range of orbits, as exists for stars in globular clus- 
ters, the inclination between the inner and outer binary or- 
bi ts is not restricte d to the coplanar case. The method given 
m IMardlind (2008) for calculating the resonance widths does 
not include the effect of the relative inclination I. The effect 
of inclination on the resonance width is under development 
jMardlindboTl and we present here a summary of the in- 
clination factors relevant to this paper. In addition to the 
inclination, the phase of the orbit also changes over time due 
to the non-Keplerian motion in the centre (see equation 1231 
for the GC potential model). For simplicity the maximum 
possible resonance width is adopted which is equivalent to 
taking the resonance angle as zero. 

From the formulation given in IMardlind J2008I) and in- 
cluding the inclination terms from IMardlind (|201 lh the res- 
onance width is given by 
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where for the leading quadrupole term in the spherical ex- 
pansion of the disturbing function one can write 
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with C22 = 3/8 and c|o = 1/4. The approximate dependence 
on the inner eccentricity is given by the functions 
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which is valid for ^ ei ^ 1. The errors between this ap- 
proximate expression and the exact integral expr ession are 
less t han 1% for a < 0.8 and < 0.1% for e ; < 0.63 (|Mardlind 
2008). The dependence of equation <(3j on the outer eccen- 
tricity is approximated by the asymptotic expression 
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As a system evolves on a secular timescale the resonance 
widths vary and are at their maximum when e, is maximum. 
Therefore the inner eccentricity in equation Q needs to be 
modified for induced eccentricity in the inner orbit due to 
the outer orbit and the secular octopole term. The maximum 
eccentricity that can be dynamically induced in the inner 
eccentricity by the outer orbit following one passage is given 
by 

eT d = [ e! (0) 2 - 2/?„ ei (0) sin(02ni) + Pi] 1/2 (11) 
where ei(0) denotes the initial inner eccentricity, and 
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where F^ 2 \e ) is given by equation ([5]). The eccentricity 
correction due to the secular octopole term, which is non- 
zero if mi ^ 7712, is 



where 
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and in the limit e; << 1 
eq (5/4)e m 3 (mi — m 2 )(ai/a ) 2 a(l 
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which is sufficiently accurate to determine the stability 
boundary for all values of e% throughout this paper. 

For inclined systems an additional secul ar effect is im- 
porta nt, this is known as the Kozai effect (see llnnanen et all 
1 1997T ) and involves a relationship between the eccentricity 
and inclination such that the maximum eccentricity induced 
by the Kozai mechanism is 
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+ 
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which all depend on the initial eccentricity of the inner bi- 
nary, ei(0), the initial relative inclination between the inner 
and outer orbits, /, and vj — vj — vji. 

The eccentricity determined by Equation [16] is the max- 
imum possible eccentricity that comes out o f this Kozai cycl e 
that gives the maximum resonance width (|Mardlindl2"01ll) . 
This means that the maximum eccentricity induced by the 
Kozai mechanism (ejf) must also be included in the inner 
eccentricity functions given in equation Q . This is achieved 
by replacing a in these equations with the theoretical max- 
imum inner eccentricity given by 
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where e\ nd is the induced eccentricity due to the outer binary 
orbit fequation lll[) and e° ct is the eccentricity correction due 
to the octopole term fequation ll3p . As the Kozai mechanism 
relates the eccentricity and inclination then the inclination 
used in equation (ff|)- (|10|) must be replaced by the maximum 
possible inclination over a Kozai cycle (Ik)- The maximum 
inclination is given by 

cos(J A -) = A (21) 



s/T 



\(Ik) = yl - cos 2 I K - 



(22) 



Each resonance width is calculated using the maximum d 
( equation 120 and inclination Ik (equations 1211 and I22[) to- 
gether with equation (|13[) . For the application to GC orbits 
in a galactic tidal field a given set of system parameters (mi, 
ma, m,3, e and Rp) are fixed and the stability for each par- 
ticular set of ei(0), a and inclination / is determined. 



3 STABILITY IN GLOBULAR CLUSTERS 

We next use the stability analysis outlined previously to 
examine how stability depends on the proximity of a star 
to the centre of the cluster. The MSC consists of a stability 
analysis to determine particular orbital configurations of the 
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Figure 1. Conceptualisation of the stability of stellar orbits in 
a globular cluster. The distance from the cluster centre associ- 
ated with the transition from unstable (dark shading) to stable 
orbits (unshaded inner region) is indicated by the ratio of outer 
to inner periods a u . The region where orbits can be found in ei- 
ther unstable or stable configurations is shown as a light shading 
between cr m i n and a m ax- The cluster is truncated at the maxi- 
mum theoretical tidal radius which is equivalent to RKing given 
by equation ((34}. 



three masses which are unstable to the escape of one of the 
masses. In the case of a star-cluster-galaxy system the only 
energetically possible end state for such an unstable config- 
uration is the escape of the star from the cluster. Therefore 
stars on orbits that are found to be unstable are predicted 
to eventually escape the cluster. 

Once again mass potentials are assumed for the galaxy 
and cluster, which is valid for the star-cluster and cluster- 
galaxy distances of interest here. The three body system is 
then described by a star rrii = 1 M©, the globular clus- 
ter M c = 10 6 M© and the galaxy M G = lO n M . In the 
notation of the stability analysis in Section the orbit of 
the star-cluster is referred to as the inner orbit composed of 
mi = Mc and m-2 = m, , and the orbit of the cluster-galaxy 
as the outer orbit (where m3 = Mg)- 

Using the MSC a boundary between predominately un- 
stable and stable orbits is sought, separated by a period ra- 
tio a u . The transition between the unstable exterior of the 
cluster and stable interior is characterised by two additional 
period ratio values <j m in and <J max . The difference between 
o~max and Omin represents the width of the transition region 
and can be used to estimate a maximum and minimum tidal 
radii range for globular clusters. A conceptualisation of the 
different stability regions for stellar orbits within a cluster 
is shown in Figure [TJ 

To determine the stability boundary inside a cluster we 
need to average over the eccentricity and inclination for the 
star-cluster orbit. In the case of the relative inclination be- 
tween the star-cluster orbit and the cluster-galaxy orbit this 
is trivial, assuming a cos I distribution. However the eccen- 
tricity distribution is not so straight forward. The approx- 
imation of such a distribution function is taken up in the 
next section. 



4 Gareth F. Kennedy 




0.2 0.4 0.6 0.8 1 



Figure 2. The distribution of eccentricities for stellar orbits in 
a Plummer sphere with 6=1 and using an isotropic velocity 
distribution. A Beta distribution is fit to the data of the form 
given in equation l|24|l and is shown as a grey curve for the best 
fit. The eccentricity distribution has a mean of = 0.47 for both 
the data and the fitting function. 



bution was sought, and the Beta distribution was found to 
provide the required fit. The Beta distribution is given by 

/(e0 = ^ y er 1 (l-e ! r i (24) 

where B(a,fJ) is the Beta function and the mean of the 
distribution is given by a /(a + 0). For the fitting function 
plotted as a grey curve in Figure [2] a = 2.691732, fj = 3, 
and B(a, /?) = 0.042898. These values ensure the same mean 
eccentricity value of e; = 0.47 for the fitting function as for 
the distribution. The cumulative probability distribution for 
equation (|24|l is 

F (ei) = „/ ^ (0.3715086? - 0.541751e? +1 
B(a,p) 

+ 0.213141e° +2 ) , (25) 

which is used when averaging over the eccentricity distri- 
bution to derive the fraction of unstable orbits in the next 
section. 



4 ECCENTRICITY DISTRIBUTION OF STARS 
IN AN ISOLATED PLUMMER SPHERE 

The MSC requires an inner eccentricity to determine the 
stability for a given system. For the orbits inside a globular 
cluster this eccentricity distribution is not obvious. In this 
section a simple simulation of N = 10 4 particles is made to 
estimate a realistic eccentricity distribution. This number 
of particles was found to give sufficient resolution for the 
required distribution, while being computationally efficient. 

A globular cluster is modelled usin g a Plummer sphere 
with gravitational potential given by (Bi nnev fc Tremainel 
1 19871 ) 
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where G is the gravitational constant, Mc is the mass of the 
cluster, r is the radial distance, and b is a parameter chosen 
to describe the compactness of the cluster. For 6 = the 
Plummer potential reduces to the potential for a point mass 
of mass Mc ■ For simplicity a value of b = 1 is assumed from 
herein. 

Numerically modelling such a cluster requires N parti- 
cles to be distributed such that the combined gravitational 
potential is equivalent to the Plummer potential, and the ve- 
locities satisfy a Maxwellian distribution. The distribution 
of particles is achieved using the mass enclosed within a par- 
ticular radius along with von Neumann's rejection technique 
for random sampling from a distribution to determine the 
veloci ties. Here the same procedure as that of lAarseth et al.l 
l|l974h is followed, except for the alterations made to account 
for the cluster compactness parameter b and for truncating 
the cluster at the King radius ( equation | 34p. Fu rther details 
of this method can be found in iKenn cdv (200l). 

To cal culate the orbit o f each particle a Bulirsch-Stoer 
integrator l|Press et al.lll986h was used until a minimum and 
maximum distance could reliably be determined, which was 
then used to give an eccentricity. The distribution of eccen- 
tricities for particle orbits in a cluster is shown in Figure [2] 

An integrable function that fits the eccentricity distri- 



5 APPLICATION TO GLOBULAR CLUSTERS 
AND APPROXIMATING THE TIDAL 
RADIUS 

To demonstrate the stability analysis process the MSC is ap- 
plied to a co-planar system with outer eccentricity e = 0.5 
following the method given in Section [2] The resonance 
widths for n:l type resonances are determined from equa- 
tion |T|) and are shown in Figure [3] (a) as a function of 
a — Vijv = T /Ti for rnilmx = 10~ 6 and rm/mi = 10 5 . 

Regions where the system resides in a single resonance 
are shaded green in Figure [3] (a) and the boundary of this 
region (the separatrix) is indicated by a black curve. The 
resonance width calculated by equation {TJ is the distance 
of the separatrix from exact resonance (n : 1). Regions where 
two or more resonances overlap are shaded in red to indicate 
theoretically unstable orbits. In the context of globular clus- 
ters, stars on these orbits are expected to eventually escape 
from the cluster. 

Unstab le systems can still occur near the separatrix (see 
Figure 15 of lMardlinell2008r ) . which means that the predicted 
unstable regions are a conservative estimate. The criterion 
of an unstable system being any system that resides in two 
resonances simultaneously - this is adopted as a quick diag- 
nostic and gives a good estimate as to where most unstable 
regions of a-ei space occur. 

By averaging over the inner eccentricity using equa- 
tion (|25[) the fraction of unstable orbits as a function of a can 
be determined. This averaging is realised by numerically de- 
termining the minimum (e m i n ) and maximum (e max ) values 
for the unstable region shown in Figure [3] (a) for a particular 
a. The fraction of unstable orbits predicted for a coplanar 
system with a fixed outer eccentricity is then 



fu n stable{o) = Pr(e 



(26) 



where F(e) is the cumulative distribution function for the 
eccentricities of particles in a Plummer sphere given by equa- 
tion (|25|) , The fraction of unstable orbits for the co- planar 
e = 0.5 system with <f> n = is shown in Figure [3] (b). The 
dashed red line in this figure indicates a u , which is the lowest 
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(a) Regions of overlapping resonances 




(b) Fraction of unstable orbits 




Figure 3. Resonance widths for a particle orbiting the cluster 
core with a perturbing galaxy on a coplanar e = 0.5 orbit for 
rri2/mi = 10 — 6 and m^/mi = 10 5 . Regions of resonance overlap 
are shaded in red to show the predicted unstable orbits in panel 
(a), while regions inside a single resonance width are shaded in 
green. Panel (b) shows the predicted fraction of unstable orbits 
as a function of the ratio of outer to inner orbital periods after 
averaging over the distribution of eccentricities in a cluster (shown 
in FigurefS}- The dashed line at a = 16 shows a u while the dotted 
lines show a min and cr max . 



a value where the fraction of unstable orbits drops beneath 
10%. This value is a representative value for characterising 
the transition from stable orbits inside the cluster to unsta- 
ble orbits in the outer regions of the cluster, as illustrated 
in Figure [TJ In Section [51 the stability boundary a u is used 
to estimate the tidal radius, after averaging over inclina- 
tion. The dotted line in Figure [3] (b) at a m in ~ 15 indicates 
the lowest a value where f uns t a ble < 0.95, while the second 
line at a max ~ 19.6 indicates the greatest a value where 
funstable > 0.05. Together these dotted lines show the width 
of the transition from stable to unstable orbits in units of 
the period ratio a for a given e (the lightly shaded region 
in Figure [TJ. 

The final fraction of unstable orbits as a function of 
the period ratio a is determined by averaging equation ()26[) 
across the range of relative inclinations between the star- 
cluster and cluster-galaxy orbits. For orbits with relative 
inclination I the inclination terms listed in Section [2] are 



used to calculate the resonance width and hence the stability 
of orbits. 

In this case one expects that cos(7) is uniformly dis- 
tributed between -1 and 1 and so the probability distribution 
function for the relative inclination can be written as 



g(I) = i(l + cos2J) 



(27) 



where ^ I ^ n. The probability of the inclination being in 
the range Ij — AI /2 ^ I ^ Ij + AI /2 is then given by 



g(l)dl 



except at the end points where it is given by 
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(29) 
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where a resolution of Ai = 7r/6 is used here. The res- 
onance widths for determining the stability of the star- 
cluster-galaxy system are now calculated using equations {JJ 
and ([2]), which includes additional terms to allow for a vari- 
able inclination. 

Using the probability distribution given by equa- 
tions (1281) to (1301). the fraction of unstable orbits becomes 



E 

3=1 



[F(e max ) - F(e min )] P(Ij) 



(31) 



where Ni nc is the total number of inclination values (set at 
12 here) spaced in increments of AI = 2n/Ni nc and F(e), 
e-max and e m ; n have been defined in the previous section. 

The fraction of unstable orbits after averaging over the 
relative inclination and eccentricities of stellar orbits within 
a cluster, as determined by equation ([31}, is shown in Fig- 
ure U (a) to (c) for e — 0.2, 0.5 and 0.8 respectively. The 
dashed vertical line in this figure again indicates the lowest 
a value where the fraction of unstable orbits drops beneath 
10% (o~ u )- The dotted lines indicate o" m i„, the lowest a value 
for which f unstable < 0.95, and a ma x, the highest a value for 

which funstable > 0.05. 

From FigureUthe transition from unstable to stable or- 
bits (a u ) increases as outer eccentricity (e ) increases. This 
is due to the dependence of the resonance width on the 
combination n£(e ) in equation (JSJ), where n refers to the 
n : 1 resonance. This quantity is always positive and as it 
increases the resonance width rapidly falls to zero. Unstable 
systems therefore require n£(e D ) to be as close to zero as 
possible, which is achieved for high values of n (and a) if e D 
is also high. Physically, this reflects the fact that an expo- 
nentially small amount of energy is exchanged between the 
inner (star-cluster) and o uter (cluster-g alaxy) orbits when 
their orbits are very wide l|Mardlindl200a '). 

The width of the transition from unstable to stable or- 
bits also increases with e as seen in the progression of panels 
(a) through (c) of Figure [4] Note that the basic structure of 
peak stability occurring at integer values of a and stability 
increasing with increasing a is consistent for all eccentric- 
ities. This phenomenon is expected since resonance widths 
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(a) e = 0.2 




Figure 4. The effect of the outer eccentricity on the fraction of 
unstable orbits against the period ratio (ex). The transition from 
unstable to stable orbits a u is shown as a vertical dashed line and 
is the lowest a value where the fraction of unstable orbits drops 
beneath 10%. The minimum and maximum cr values associated 
with the width of this transition are shown as dotted lines. Note 
the increase in the range of a as e D increases. 

from the n : 1 and n + 1 : 1 resonances (where n < a < n + 1) 
overlap at the midpoint between these resonances, as seen 
previously in Figure [3] (a) for the coplanar case. 

The results for all eccentricity values are shown in Fig- 
ure[Sl from which a near exponential dependence of a u on e D 
is seen. This is also true for the width of the transition be- 
tween unstable and stable orbits as represented by a min / max 
and shown as dotted curves on each side of a u in Figure [5] 
Note that the coarseness in the curves for all a values for 
low eccentricities is due to a resolution of 1 cr used when 
producing this figure. For high eccentricities it is clear from 
the figure that there is a large range in period ratios from 
which stars can potentially escape the cluster. 

The transition value of a u will be used in the next sec- 
tion to calculate the tidal radius for a given perigalacticon 
and eccentricity of the cluster orbit about the galaxy. The 
width of the transition from unstable to stable orbits, given 
by &min/max, wm be used to provide approximate error bars 
associated with the tidal radius. 

This tidal radius calculation is intended to be generally 
applicable to the entire system of globular clusters. However 
in determining <j(e) the mass ratios have been taken as con- 
stant. This is implicitly reflected in the stellar cluster model 
used to determine the probability distribution for the eccen- 
tricity of stellar orbits within the cluster (equation (|25[) ), 
This cluster consisted of IMq particles in a cluster of mass 
Mc = 1O 6 M0, which will not be true of all clusters. 

Changing the cluster mass is expected to have a neg- 
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e out 

Figure 5. The effect of the galactic orbital eccentricity of the 
cluster on the transition from unstable to stable star-cluster or- 
bits. The transition a u is shown as a solid line and is used as 
an estimate for the tidal radius, while the associated transition 
width values cr min / max are shown as dotted lines. 

ligible effect on the resonance width calculation and there- 
fore on the stability boundary as a function of eccentricity. 
This can be demonstrated by considering the mass depen- 
dence of equation ([3]) when rri2 = nij <C Mc = m\ and 
mi = Mc <H Mc = ra-i, which can be simplified to 

where for distant globular clusters a ~ 10 kpc, aj ~ 10 pc 
and Mc = lO n M means that this term is independent of 
mass for 10 4 < Mc/Mq < 10 6 . We can therefore apply the 
cr(e) relationship derived in this section to all cluster masses 
for distant globular clusters. 



6 COMPARISON WITH PREVIOUS 
THEORETICAL WORK 

Two theoretical estimates for the tidal radius from the liter- 
ature are compared to the tidal radius derived from the MSC 
method. The first estimate from the literature is the most 
commonly u sed tidal radius estimate in the field derived by 
iKind (|l962h and will be referred to as the King radius. The 
second is an extended analytical determination that was also 
compared to N-body simulations for a s ingle set of orbita l 
parameters, this determination is given in I Read et al.l (2006) 
and will be referred to as the Read radius. 

To aid comparison between these estimates the eccen- 
tricity dependence is separated out so that the tidal radius 
can be written as 

rt =MII) 1/3/(e) (33) 

where Mc and Mc are the masses of the cluster and galaxy 
respectively, e is the eccentricity of the clusters orbit around 
the galaxy, and R p is the distance of closest approach to the 
galaxy, referred to as the perigalacticon. 

The simplest case to determine analytically is to con- 
sider a star located where the acceleration on the star in the 
rotating frame is zero and the velocity of the star relative to 
the cluster centre is also zero. Such a star will be on a radial 
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orbit with respect to the centre of the cluster. For a star on 
a radial orbit and using point mass potentials for the cluster 
and galax y the eccen tricity dependence of the tidal radius is 
given by (|Kinglll962r ) 

/(e) = fc(3 + e)- 1/3 (34) 

where the constant k ~ 0.7 was introduced by iKeenanl 
(1981) to better fit observations of the galactic globular clus- 
ters. 

Recently the analysis o f iKind (|l962h has been extended 
to star s on circular orbi ts on either prograde or retrograde 
orbits (|Read et al J 12006). A star has zero acceleration in a 
coordinate frame rotating with the motion of a point mass 
cluster if the eccentrici ty dependence of the distance from 
the cluster is given by (|Read et al.l l2006) 



1 



1/3 




+ 1 + 



2/3 



1 + — 

^ 1+e 



(35) 



where a — denotes a star on a radial orbit inside the clus- 
ter and reduces to equation (|34l) in this case. Non-radial mo- 
tion is restricted to the case of stars on circular orbits and 
is described in the tidal radius equation by setting a = 1 
or —1 for prograde and retrograde orbits respectively (note 
that a = — 1 gives r t greater than the King radius). For 
later comparison we define the Read radius as equation (|35p 
with a = 1, representing an easily tidally stripped cluster. 
The Read radius compared well with two N-body simula- 
tions of 1O 7 M0 satellite clusters using N = 10 5 particles 
with orbital parameters of perigalacticon Rp/Rx/2 = 267 
and eccentricity e = 0.0 an d perigalacticon R p/R 1 / 2 = 77 
and eccentricity e = 0.57 |Read et al.ll2006t ). A summary 
of alternate equations for t he tidal radius d esigned to fit 
observations can be found in iBellazzinil (|2004l '). 

For the stability boundary determined in Section [2] to 
be useful it must first be converted into an equivalent radius 
from the cluster centre. Assuming that the gravitational po- 
tential of the cluster is well approximated by a point mass 
at distances of the tidal radius for clusters of interest then 
the tidal radius can be written in the same form of equa- 
tion (|33[) . By converting the period ratio shown in Figure [5] 
into an equivalent semi- major axis (oi) via T = oTi and 
using the resulting a.i as r t then the eccentricity dependence 
itself as defined by /(e) in equation (|33[) can be determined. 
The values for /(e) as determined from the period ratio a 
dependence on eccentricity (Figure [S| are shown as the data 
points in Figure [6] These data can be fit by 



fie) 



exp 



N=7 

E< 

i=0 



(36) 



where the coefficients are given in Table [T] for the minimum, 
maximum and indicative (chaos) radii. 

From Figure [5] we see that the boundary between sta- 
ble and unstable orbits occurs interior to both the Read and 
King tidal radius estimates. However the stability boundary 
is not perfectly equivalent to the tidal radius. This is because 
although chaotic orbits will eventually result in the escape of 
the star from the globular cluster, this occurs via a random 
walk process. Therefore there will be stars remaining out- 
side the distance determined to be the stability boundary, 
sometimes for many crossing times. These stars will have a 





Rchaos 


Rhlin 


RMax 


on 


-0.897753 


-0.892899 


-0.888868 


"1 


-2.28272 


-1.83631 


-1.93401 


«2 


16.8616 


10.3382 


13.844 


«3 


-74.0133 


-45.1384 


-63.0546 


(14 


177.136 


108.206 


163.39 


a 5 


-236.158 


-144.186 


-239.155 


do 


163.999 


99.5267 


182.311 


a 7 


-46.3734 


-27.6654 


-56.3393 



Table 1. Coefficients for fit to Rchaos and the minimum and 
maximum extents of the marginally chaotic zone. For radii r < 
Rutin a U orbits are expected to be stable, whereas for r > Rmux 
they are expected to be chaotic. 




Figure 6. The eccentricity dependence of the tidal radii asso- 
ciated with the MSC prediction are shown as data points along 
with a solid line fit given by equation H36I I. the King radius as 
dashed lines and the Read radius as dotted lines. Coefficients to 
fitting functions given in Table [T] along with the minimum (red) 
and maximum (green) fits. 



different velocity dispersion profile to what is expected for a 
relaxed and isolated globular cluster, a fact that forms the 
basis for the second paper in this series (|Kennedvll26TTh . 

The dangers of assuming that the tidal radius acts as 
an instant remover of sta rs has also been pointed out by 
iFukushige fc Heggid (|2000h . They found for GCs on circular 
orbits that the escape timescales for stars beyond the tidal 
radii could be long enough to allow some stars to stay in this 
region indefinitely. This means that all predicted tidal radii 
will be lower limits, as stars outside this region can still be 
close to the GC while remaining formally unbound. 



7 DISCUSSION AND CONCLUSIONS 

The tidal radius determined here differs from previous esti- 
mates in that it emerges naturally from a stability analysis 
of the general three-body problem. It was found that the 
eccentricity of the cluster-galaxy orbit has a far more sig- 
nificant effect on the tidal radius of globular clusters than 
previous theoretical results have found. 

The approach adopted here means that what is actu- 
ally predicted is the boundary between stable and unstable 
orbits for stars inside a cluster, which is interior to previous 
values for all parameters of the cluster-galaxy orbit. This 
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was found to be the case i n com parison to the commonly 
used tidal radius by King ( 19621 ) and to the more recent 



radius by iRead etal 



|2ooa T One key difference between 
the tidal radius estimation presented here and previously 
is that a much stronger dependence on the eccentricity of 
the cluster-galaxy orbit is predicted here. This dependence 
can be studied using N-body simulations, which would also 
give insight into how long stars on unstable orbits take to 
random walk out of the cluster. 

From a practical point of view a major outcome of this 
work is the derivation of an easy to use tidal radius of the 
form 

r *=*{mT™ (37) 

where /(e) is given by a function (equation I36|) fitted to 
stability results determined analytically from the stability 
of the general three-body problem. This has a wide range of 
applications, including the effect of unstable orbits on the 
velocity dispersion p rofile for Milky Way globular clusters, 
which is taken up in iKenne 



8 ACKNOWLEDGEMENTS 

I am indebted to the three-body stability work of Rosemary 
Mardling, in particular for supplying the unpublished incli- 
nation terms used in this study. 



REFERENCES 

Aarseth S. J., Henon M., Wielen R., 1974, A&A, 37, 183 

Bellazzini M., 2004, MNRAS, 347, 119 

Binney J., Tremaine S., 1987, Galactic dynamics. Prince- 
ton, NJ, Princeton University Press, 1987 

Fukushige T., Heggie D. C, 2000, MNRAS, 318, 753 

Innanen K. A., Zheng J. Q., Mikkola S., Valtonen M. J., 
1997, AJ, 113, 1915 

Keenan D. W., 1981, A&A, 95, 340 

Kennedy G., 2008, Ph.d. thesis, Monash University 

Kennedy G. F., , 2011, this volume 

King I., 1962, AJ, 67, 471 

Mardling R. A., 2008, in Aarseth S. J., Tout C. A., 
Mardling R. A., eds, Lecture Notes in Physics, Vol. 760: 
The Cambridge N-body Lectures Resonance, chaos and 
stability: the three-body problem in astrophysics 

Mardling R. A., , 2011, in preparation 

Press W. B., Flannery B. P., Teukolsky S. A., Vettering 
W. T., 1986, Numerical Recipes: The Art of Scientific 
Computing. Cambridge Univ. Press, Cambridge 

Read J. I., Wilkinson M. 1., Evans N. W., Gilmore G., 
Kleyna J. T., 2006, MNRAS, 366, 429 



